Imports MatrixDataStructures

Public Class BasicPixelSizeMatrixInterpolation
    Public smooth As Gaussian3D
    Public Function ResizePixels(ByVal OriginalMatrix As Matrix3DInt16, ByVal FinalX As Double, ByVal FinalY As Double, ByVal FinalZ As Double) As Matrix3DInt16

        Dim TempResult As New Matrix3DInt16(FinalX, FinalY, FinalZ)

        'suppose final is 10 10 10
        'original is 5 5 5
        'then ratio is 2 2 2
        'so when looking for 10 10 10 value, you are lookign for 5 5 5 value
        'ie value for i/ratiox j / ratioy k / raioz value
        Dim RatioX As Double = FinalX / OriginalMatrix.x
        Dim RatioY As Double = FinalY / OriginalMatrix.y
        Dim RatioZ As Double = FinalZ / OriginalMatrix.z

        Dim MaxOx, MaxOy, MaxOz As Integer
        MaxOx = OriginalMatrix.x - 1
        MaxOy = OriginalMatrix.y - 1
        MaxOz = OriginalMatrix.z - 1

        Dim temp1, temp2, temp3 As Integer
        For i As Integer = 0 To FinalX - 1
            For j As Integer = 0 To FinalY - 1
                For k As Integer = 0 To FinalZ - 1

                    'check if the ratio isnt out of the original matrix
                    temp1 = CInt(IIf(i / RatioX < MaxOx, i / RatioX, MaxOx))
                    temp2 = CInt(IIf(j / RatioY < MaxOy, j / RatioY, MaxOy))
                    temp3 = CInt(IIf(k / RatioZ < MaxOz, k / RatioZ, MaxOz))

                    TempResult(i, j, k) = OriginalMatrix(temp1, temp2, temp3)
                Next
            Next
        Next

        Return TempResult
    End Function
    Public Function ResizePixelsReturnSingle(ByVal OriginalMatrix As Matrix3DInt16, ByVal FinalX As Double, ByVal FinalY As Double, ByVal FinalZ As Double) As Matrix3DSingle

        Dim TempResult As New Matrix3DSingle(FinalX, FinalY, FinalZ)

        'suppose final is 10 10 10
        'original is 5 5 5
        'then ratio is 2 2 2
        'so when looking for 10 10 10 value, you are lookign for 5 5 5 value
        'ie value for i/ratiox j / ratioy k / raioz value
        Dim RatioX As Double = FinalX / OriginalMatrix.x
        Dim RatioY As Double = FinalY / OriginalMatrix.y
        Dim RatioZ As Double = FinalZ / OriginalMatrix.z

        Dim MaxOx, MaxOy, MaxOz As Integer
        MaxOx = OriginalMatrix.x - 1
        MaxOy = OriginalMatrix.y - 1
        MaxOz = OriginalMatrix.z - 1

        Dim temp1, temp2, temp3 As Integer
        For i As Integer = 0 To FinalX - 1
            For j As Integer = 0 To FinalY - 1
                For k As Integer = 0 To FinalZ - 1

                    'check if the ratio isnt out of the original matrix
                    temp1 = CInt(IIf(i / RatioX < MaxOx, i / RatioX, MaxOx))
                    temp2 = CInt(IIf(j / RatioY < MaxOy, j / RatioY, MaxOy))
                    temp3 = CInt(IIf(k / RatioZ < MaxOz, k / RatioZ, MaxOz))

                    TempResult(i, j, k) = OriginalMatrix(temp1, temp2, temp3)
                Next
            Next
        Next

        Return TempResult
    End Function
    'resize pixels is done to resize the matrix to its OWN fov
    Public Function ResizePixels(ByVal OriginalMatrix As Matrix3DSingle, ByVal FinalX As Integer, ByVal FinalY As Integer, ByVal FinalZ As Integer) As Matrix3DSingle

        Dim TempResult As New Matrix3DSingle(FinalX, FinalY, FinalZ)

        'suppose final is 10 10 10
        'original is 5 5 5
        'then ratio is 2 2 2
        'so when looking for 10 10 10 value, you are lookign for 5 5 5 value
        'ie value for i/ratiox j / ratioy k / raioz value
        Dim RatioX As Double = FinalX / OriginalMatrix.x
        Dim RatioY As Double = FinalY / OriginalMatrix.y
        Dim RatioZ As Double = FinalZ / OriginalMatrix.z

        Dim MaxOx, MaxOy, MaxOz As Integer
        MaxOx = OriginalMatrix.x - 1
        MaxOy = OriginalMatrix.y - 1
        MaxOz = OriginalMatrix.z - 1

        Dim temp1, temp2, temp3 As Integer
        For i As Integer = 0 To FinalX - 1
            For j As Integer = 0 To FinalY - 1
                For k As Integer = 0 To FinalZ - 1

                    'check if the ratio isnt out of the original matrix
                    temp1 = CInt(IIf(i / RatioX < MaxOx, i / RatioX, MaxOx))
                    temp2 = CInt(IIf(j / RatioY < MaxOy, j / RatioY, MaxOy))
                    temp3 = CInt(IIf(k / RatioZ < MaxOz, k / RatioZ, MaxOz))

                    TempResult(i, j, k) = OriginalMatrix(temp1, temp2, temp3)
                Next
            Next
        Next
        Return TempResult
    End Function
    Public Function ResizePixels(ByVal OriginalMatrix As Matrix4DInt16, ByVal FinalX As Integer, ByVal FinalY As Integer, ByVal FinalZ As Integer) As Matrix4DInt16

        Dim TempResult As New Matrix4DInt16(FinalX, FinalY, FinalZ, OriginalMatrix.l)

        'WORKS ONLY IF L COORDINATE is VOLUME/TIMEPOINT 

        'suppose final is 10 10 10
        'original is 5 5 5
        'then ratio is 2 2 2
        'so when looking for 10 10 10 value, you are lookign for 5 5 5 value
        'ie value for i/ratiox j / ratioy k / raioz value
        Dim RatioX As Double = FinalX / OriginalMatrix.i
        Dim RatioY As Double = FinalY / OriginalMatrix.j
        Dim RatioZ As Double = FinalZ / OriginalMatrix.k

        Dim MaxOx, MaxOy, MaxOz As Integer
        MaxOx = OriginalMatrix.i - 1
        MaxOy = OriginalMatrix.j - 1
        MaxOz = OriginalMatrix.k - 1

        Dim temp1, temp2, temp3 As Integer
        For l As Integer = 0 To OriginalMatrix.l - 1
            For i As Integer = 0 To FinalX - 1
                For j As Integer = 0 To FinalY - 1
                    For k As Integer = 0 To FinalZ - 1

                        'check if the ratio isnt out of the original matrix
                        temp1 = CInt(IIf(i / RatioX < MaxOx, i / RatioX, MaxOx))
                        temp2 = CInt(IIf(j / RatioY < MaxOy, j / RatioY, MaxOy))
                        temp3 = CInt(IIf(k / RatioZ < MaxOz, k / RatioZ, MaxOz))

                        TempResult.data(i, j, k, l) = OriginalMatrix.data(temp1, temp2, temp3, l)
                    Next
                Next
            Next
        Next
        Return TempResult
    End Function

    'resize pixels is done for matching
    Public Function ResizeVolumes(ByVal myMatrix As Matrix3DInt16, ByVal myCenter As CoordinateD, ByVal targetCenter As CoordinateD, ByVal myFOV As CoordinateD, ByVal targetFOV As CoordinateD) As Matrix3DInt16

        Dim returnMatrix As New Matrix3DInt16(Math.Round(targetFOV.x), Math.Round(targetFOV.y), Math.Round(targetFOV.z))
        Dim matrixc, absolutec As New CoordinateD(0, 0, 0)
        For i As Integer = 0 To myMatrix.x - 1
            For j As Integer = 0 To myMatrix.y - 1
                For k As Integer = 0 To myMatrix.z - 1
                    'calculate the absolute cooridnate for (i,j,k) in actual 3d space
                    absolutec = ToAbsoluteCoordinate(New CoordinateD(i, j, k), myFOV, myCenter)
                    'find its indexes in the virtual 3d space
                    matrixc = ToMatrixCoordinate(absolutec, targetFOV, targetCenter)
                    If matrixc.x > 0 And matrixc.y > 0 And matrixc.z > 0 And _
                        Math.Round(matrixc.x) < returnMatrix.x And Math.Round(matrixc.y) < returnMatrix.y And Math.Round(matrixc.z) < returnMatrix.z Then

                        returnMatrix(Math.Round(matrixc.x), Math.Round(matrixc.y), Math.Round(matrixc.z)) = myMatrix(i, j, k)
                    End If
                Next
            Next
        Next
        Return returnMatrix
    End Function
    Public Function ResizeVolumes(ByVal myMatrix As Matrix3DSingle, ByVal myCenter As CoordinateD, ByVal targetCenter As CoordinateD, ByVal myFOV As CoordinateD, ByVal targetFOV As CoordinateD) As Matrix3DSingle
        Dim returnMatrix As New Matrix3DSingle(Math.Round(targetFOV.x), Math.Round(targetFOV.y), Math.Round(targetFOV.z))
        Dim matrixc, absolutec As New CoordinateD(0, 0, 0)
        For i As Integer = 0 To myMatrix.x - 1
            For j As Integer = 0 To myMatrix.y - 1
                For k As Integer = 0 To myMatrix.z - 1
                    'calculate the absolute cooridnate for (i,j,k) in actual 3d space
                    'MsgBox("ijk: " & i & "," & j & "," & k)
                    absolutec = ToAbsoluteCoordinate(New CoordinateD(i, j, k), myFOV, myCenter)
                    'MsgBox("abs: " & absolutec.x & "," & absolutec.y & "," & absolutec.z)
                    'find its indexes in the virtual 3d space
                    matrixc = ToMatrixCoordinate(absolutec, targetFOV, targetCenter)
                    'MsgBox("i'j'k': " & matrixc.x & "," & matrixc.y & "," & matrixc.z)
                    If matrixc.x > 0 And matrixc.y > 0 And matrixc.z > 0 And _
                        matrixc.x < returnMatrix.x And matrixc.y < returnMatrix.y And matrixc.z < returnMatrix.z Then
                        returnMatrix(Math.Round(matrixc.x), Math.Round(matrixc.y), Math.Round(matrixc.z)) = myMatrix(i, j, k)
                    End If
                Next
            Next
        Next
        Return returnMatrix
    End Function

    Private Function ToAbsoluteCoordinate(ByVal myMatrixCoordinate As CoordinateD, ByVal myfov As CoordinateD, ByVal myrealcenter As CoordinateD) As CoordinateD
        Return New CoordinateD(myrealcenter.x - (myfov.x / 2) + myMatrixCoordinate.x, _
myrealcenter.y - (myfov.y / 2) + myMatrixCoordinate.y, _
myrealcenter.z - (myfov.z / 2) + myMatrixCoordinate.z)
    End Function
    Private Function ToMatrixCoordinate(ByVal myAbsoluteCoordinate As CoordinateD, ByVal myfov As CoordinateD, ByVal myrealcenter As CoordinateD) As CoordinateD
        Return New CoordinateD(-1 * myrealcenter.x + (myfov.x / 2) + myAbsoluteCoordinate.x, _
-1 * myrealcenter.y + (myfov.y / 2) + myAbsoluteCoordinate.y, _
-1 * myrealcenter.z + (myfov.z / 2) + myAbsoluteCoordinate.z)
    End Function

End Class

'copied from segmentation routine
Public Class Gaussian3D
    Dim Variance As Double
    Dim HalfKernelSize As Integer

    Sub New()
    End Sub
    Sub New(ByVal KernelSize As Integer, ByVal Variance As Double)
        If KernelSize Mod 2 = 1 Then
            Me.HalfKernelSize = CInt((KernelSize - 1) / 2)
        Else
            Me.HalfKernelSize = CInt(KernelSize / 2)
        End If
        Me.Variance = Variance
    End Sub

    Private Function CreateKernel() As Matrix3DSingle
        Dim Kernel As New Matrix3DSingle(HalfKernelSize * 2 + 1, HalfKernelSize * 2 + 1, HalfKernelSize * 2 + 1)
        Dim val As Double = 0.0

        For i As Integer = -1 * HalfKernelSize To 1 * HalfKernelSize
            For j As Integer = -1 * HalfKernelSize To 1 * HalfKernelSize
                For k As Integer = -1 * HalfKernelSize To 1 * HalfKernelSize
                    val = (i * i + j * j + k * k) / (-2 * Variance)
                    Kernel(i + HalfKernelSize, j + HalfKernelSize, k + HalfKernelSize) = CSng(Math.Exp(val) / ((2 * Math.PI * Variance) ^ 0.5))
                Next
            Next
        Next
        Return Kernel
    End Function

    Public Function ApplyGaussian(ByVal DataMatrix As Matrix3DInt16) As Matrix3DInt16
        Dim Kernel As Matrix3DSingle = CreateKernel()
        Dim temp As New Matrix3DInt16(DataMatrix.x, DataMatrix.y, DataMatrix.z)
        For i As Integer = 0 To DataMatrix.x - 1
            For j As Integer = 0 To DataMatrix.y - 1
                For k As Integer = 0 To DataMatrix.z - 1
                    If DataMatrix.data(i, j, k) > 0 Then
                        'now for each point in the datamatrix
                        'start fresh from the left-topmost of the kernel
                        For ki As Integer = -1 * HalfKernelSize To 1 * HalfKernelSize
                            For kj As Integer = -1 * HalfKernelSize To 1 * HalfKernelSize
                                For kk As Integer = -1 * HalfKernelSize To 1 * HalfKernelSize
                                    'skip points that are outside the image when the kernel convolves 
                                    'with the data matrix
                                    If i + ki < 0 Or i + ki > DataMatrix.x - 1 Or _
                                    j + kj < 0 Or j + kj > DataMatrix.y - 1 Or _
                                    k + kk < 0 Or k + kk > DataMatrix.z - 1 Then
                                        'tempvalue = 0
                                        Exit For
                                    End If
                                    temp(i, j, k) += DataMatrix.data(i + ki, j + kj, k + kk) * Kernel.data(ki + HalfKernelSize, kj + HalfKernelSize, kk + HalfKernelSize)
                                Next
                            Next
                        Next
                    End If
                Next
            Next
        Next
        Return temp
    End Function
    Public Function ApplyGaussian(ByRef DataMatrix As Matrix3DSingle) As Matrix3DSingle
        Dim Kernel As Matrix3DSingle = CreateKernel()
        Dim temp As New Matrix3DSingle(DataMatrix.x, DataMatrix.y, DataMatrix.z)

        For i As Integer = 0 To DataMatrix.x - 1
            For j As Integer = 0 To DataMatrix.y - 1
                For k As Integer = 0 To DataMatrix.z - 1
                    'now for each point in the datamatrix
                    If DataMatrix(i, j, k) > 0 Then
                        'start fresh from the left-topmost of the kernel
                        For ki As Integer = -1 * HalfKernelSize To 1 * HalfKernelSize
                            For kj As Integer = -1 * HalfKernelSize To 1 * HalfKernelSize
                                For kk As Integer = -1 * HalfKernelSize To 1 * HalfKernelSize
                                    'skip points that are outside the image when the kernel convolves 
                                    'with the data matrix
                                    If i + ki < 0 Or i + ki > DataMatrix.x - 1 Or _
                                    j + kj < 0 Or j + kj > DataMatrix.y - 1 Or _
                                    k + kk < 0 Or k + kk > DataMatrix.z - 1 Then
                                        'tempvalue = 0
                                        Exit For
                                    End If
                                    temp(i, j, k) += DataMatrix(i + ki, j + kj, k + kk) * Kernel(ki + HalfKernelSize, kj + HalfKernelSize, kk + HalfKernelSize)
                                Next
                            Next
                        Next
                    End If
                Next
            Next
        Next
        Return temp
    End Function
End Class